# GNU General Public License
#
# Program Tops - a stack-based computing environment
# Copyright (C) 1999-2005  Dale R. Williamson
#
# Author and copyright holder of op4_master.c, loadop4.c, saveop4.c, op4.c:
# Al Danial <al.danial@gmail.com>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#  This file is automatically generated from the tops file
#  admin/templates/op4_master.c.  Changes applied to this
#  file will be lost.
#


#include "sparse.h"
#include "op4.h"

#include "mex.h"
#define  MAX(a,b)              ((a) >  (b)) ? (a)  : (b)

def get_file_suffix(file,ext):
    """
    constants
    would like to have the constant arrays below defined in op4.h
    but they cause 'multiple definition' errors
    """

    op4_store_str = [ ['d',  'n',  0],
                      ['s',  '1',  0],
                      ['s',  '2',  0]
                    ];
    op4_words_per_term = [-1, 1, 2, 2, 4];
    Nm_per_term  = [1,  # numbers per term[real data type]    = 1
                    2]; # numbers per term[complex data type] = 2
    op4_type_str = [ ['R',  'S',  0],
                     ['R',  'D',  0],
                     ['C',  'S',  0],
                     ['C',  'D',  0]
                   ];

# internal functions

def op4_filetype(filename):
    """
    Returns
          1 if a text file
          2 if a binary file native  endian
          3 if a binary file flipped endian
          0 if none of the above or file read error
    """
    DEBUG = 0
    int   word, type, nRead
    #char *w

    #w  = (char *) &word;
    fp = open(filename, "rb"); # open in binary mode

    nRead = fread(&word, BYTES_PER_WORD, 1, fp);
    if DEBUG:
        print "op4_filetype first word is %d  chars:[%c][%c][%c][%c]\n" %(word, w[0], w[1], w[2], w[3]);
    if(nRead == 1): # in units of words
        if(word >= 538976288):
            # If this is a text file the smallest possible value for
            #  the integer comprising the first four bytes is
            #  538976288 which corresponds to four blank spaces.

            if((ISBLANK(      w[0]) and ISBLANK(      w[1]) and   # \s\s\s\s
                ISBLANK(      w[2]) and ISBLANK(      w[3])) or
               (ISBLANK(      w[0]) and ISBLANK(      w[1]) and   # \s\s\s\d
                ISBLANK(      w[2]) and isdigit((int) w[3])) or
               (ISBLANK(      w[0]) and ISBLANK(      w[1]) and   # \s\s\d\d
                isdigit((int) w[2]) and isdigit((int) w[3])) or
               (ISBLANK(      w[0]) and isdigit((int) w[1]) and   # \s\d\d\d
                isdigit((int) w[2]) and isdigit((int) w[3])) or
               (isdigit((int) w[0]) and isdigit((int) w[1]) and   # \d\d\d\d
                isdigit((int) w[2]) and isdigit((int) w[3]))):
                type = 1;
            else:
                type = 0;

        elif(word == 24):        # Fortran binary record header,
            type = 2;            # native endian, for 6 word block

        elif(word == 402653184): # Fortran binary record header,
            type = 3;            # opposite endian, for 6 word block
        else:     # the first word indicates the file is not an .op4
            type = 0;
    else:        # unable to read four bytes from the file
        type = 0;
    fp.close();
    return type;

def op4_scan(  filename,  # in
                n_mat,    # out number of matrices
                name,     # out matrix names
                storage,  # out 0=dn; 1=sp1; 2=sp2
                nRow,     # out number of rows
                nCol,     # out number of columns
                nStr,     # out number of strings
                nNnz,     # out number of nonzero terms
                type,     # out 1=RS; 2=RD; 3=CS; 4=CD
                form,     # out matrix form 6=symm, etc
                digits,   # out size of mantissa
                offset    # out byte offset to matrix
                ):
    """
    Count the number and determine storage type of matrices
    in the text op4 file by scanning for matrix, column & string headers.
    Driver to text and binary versions.
    """
    #int i, filetype, endian,result # = 0  causes segfault

    filetype = op4_filetype(filename); # 1 if a text file
                                       # 2 if a binary file native  endian
                                       # 3 if a binary file flipped endian
                                       # 0 if none of the above or file error

    if filetype==1:
        result = op4_scan_t(filename  ,  # in
                   n_mat     ,  # out number of matrices
                   name      ,  # out matrix names
                   storage   ,  # out 0=dn; 1=sp1; 2=sp2
                   nRow      ,  # out number of rows
                   nCol      ,  # out number of columns
                   nStr      ,  # out number of strings
                   nNnz      ,  # out number of nonzero terms
                   type      ,  # out 1=RS; 2=RD; 3=CS; 4=CD
                   form      ,  # out matrix form 6=symm, etc
                   digits    ,  # out size of mantissa
                   offset);     # out byte offset to matrix

    elif filetype==2 or filetype==3:
            endian = filetype - 2;
            result = op4_scan_b(filename  ,  # in
                       endian    ,  # in  0=native   1=flipped
                       n_mat     ,  # out number of matrices
                       name      ,  # out matrix names
                       storage   ,  # out 0=dn; 1=sp1; 2=sp2
                       nRow      ,  # out number of rows
                       nCol      ,  # out number of columns
                       nStr      ,  # out number of strings
                       nNnz      ,  # out number of nonzero terms
                       type      ,  # out 1=RS; 2=RD; 3=CS; 4=CD
                       form      ,  # out matrix form 6=symm, etc
                       offset);     # out byte offset to matrix
            if endian:
                for i in range(n_mat):
                    digits[i] = -1;
            else:
                for i in range(n_mat):
                    digits[i] =  0;
    else:
        result = 0;
    return result;

def  op4_scan_t(filename  ,  # in
                n_mat     ,  # out number of matrices
                name,     # out matrix names
                storage,  # out 0=dn; 1=sp1; 2=sp2
                nRow,     # out number of rows
                nCol,     # out number of columns
                nStr,     # out number of strings
                nNnz,     # out number of nonzero terms
                type,     # out 1=RS; 2=RD; 3=CS; 4=CD
                form,     # out matrix form 6=symm, etc
                digits,   # out size of mantissa
                offset    # out byte offset to matrix
               ):
    """
       Count the number and determine storage type of matrices
       in the text op4 file by scanning for matrix, column & string headers.
       Text version.
       Keep track of byte location, matrix name and matrix dims.
       Here are some typical headers:
      10      10       6       2EYE     1P,3E23.16
      10     -10       6       2EYE     1P,3E23.16
123456789x123456789x123456789x123456789x123456789x123456789x
                                 |-> optional from here on

       The first is for dense and sparse form 1, and the
       second is for sparse form 2.
       The headers are distinguished from numeric data in that
       there is no decimal point in the third column--a quality
       all text based numeric data has, as in
-3.937729118E+02-4.689545312E+07 3.766712501E+07-6.581960491E+06 0.000000000E+00
     """
    DEBUG = 0
    #char  line[OP4_TXT_LINE_SIZE + 1];
    #int   n_text_cols, text_col_width, mixed, nwords, start_row, nRead;
    #int  location;
    #char  junk1[12], junk2[12], junk3[12];

    if DEBUG:
        print "op4_scan_t file=[%s]\n" %(filename)
    n_mat = 0;
    fp = open(filename, "r");

    if DEBUG:
        print "op4_scan_t opened ok\n"
    while (not feof(fp)):
        location = fp.tell();
        fgets(line, OP4_TXT_LINE_SIZE, fp);
        if(strlen(line) >= 8):
            # actually the op4 file is illegal if line length < 8
            
            lineType = op4_line_type(line)
                if lineType==OP4_TEXT_NUMERIC:
                    break
                elif lineType==OP4_TEXT_HEADER: # {{{2
                    if(n_mat and not storage[n_mat-1]):
                        # previous matrix stored in dense scheme
                        nNnz[n_mat-1] = nRow[n_mat-1]nCol[n_mat-1];

                    storage[n_mat] = 0;  # dense unless proven otherwise
                    nStr[n_mat]    = 0;
                    nNnz[n_mat]    = 0;
                    if(line[41] == 'P'):
                        nRead = sscanf(line,
                                      "%8d%8d%8d%8d%8s%[ 1P,]%1d%1[E]%d%1[.]%d",
                                      &nCol[n_mat],
                                      &nRow[n_mat],
                                      &form[n_mat],
                                      &type[n_mat],
                                       name[n_mat],
                                       junk1,           #   1P,
                                      &n_text_cols,
                                       junk2,           #   E
                                      &text_col_width,
                                       junk3,           #   .
                                      &digits[n_mat]);
                        if(nRead != 11):
                            fp.close()
                            return 0;
                    else:
                        n_text_cols    =  5
                        text_col_width = 15
                        digits[n_mat]  =  9
                        nRead = sscanf(line, "%8d%8d%8d%8d%8s",
                                      &nCol[n_mat],
                                      &nRow[n_mat],
                                      &form[n_mat],
                                      &type[n_mat],
                                       name[n_mat])
                        if(nRead != 5):
                            fp.close()
                            return 0;

                    nRow[n_mat]   = fabs(nRow[n_mat])
                    offset[n_mat] = location
                    ++(n_mat);
                elif lineType==OP4_TEXT_NEW_COLUMN: # {{{2
                    if DEBUG:
                        print "new col: %s" %(line)
                elif lineType==OP4_TEXT_STRING_2: # {{{2
                    sscanf(line, "%8d%8d", &nwords, &start_row);
                    --nwords;
                    nNnz[n_mat-1] += nwords/op4_words_per_term[type[n_mat-1]]
                    if DEBUG:
                        print "s2:  nwords=%d  nNnz=%d\n" %(nwords, nNnz[n_mat-1]);
                    storage[n_mat-1] = 2;
                    ++nStr[ n_mat-1];
                elif lineType==OP4_TEXT_STRING_1:
                    if DEBUG:
                        print "s1:  %s" %(line)
                    sscanf(line, "%8d", &mixed);
                    nwords    = (mixed/65536) - 1;
                    start_row = mixed - 65536 * (nwords + 1);
                    nNnz[   n_mat-1] += nwords/op4_words_per_term[type[n_mat-1]];
                    storage[n_mat-1] = 1;
                    ++nStr[n_mat-1];
                elif lineType==OP4_TEXT_ERROR:
                    fp.close()
                    return 0;

    offset[n_mat] = fp.tell();  # last byte in the file
    fp.close()
    if not storage[n_mat-1]:
        # last matrix stored in dense scheme
        nNnz[n_mat-1] = nRow[n_mat-1]nCol[n_mat-1];
    return 1;


def op4_scan_b(filename, # in
                endian ,  # in  0=native   1=flipped
                n_mat ,   # out number of matrices
                name,     # out matrix names
                storage,  # out 0=dn; 1=sp1; 2=sp2
                nRow   ,  # out number of rows
                nCol   ,  # out number of columns
                nStr   ,  # out number of strings
                nNnz   ,  # out number of nonzero terms
                type   ,  # out 1=RS; 2=RD; 3=CS; 4=CD
                form   ,  # out matrix form 6=symm, etc
                offset    # out byte offset to matrix
               ):
    """
    Count the number and determine storage type of matrices
    in the text op4 file by scanning for matrix, column & string headers.
    Binary version.
    """
    DEBUG = 0
    #int   record_length, is_header, result  # = 0  causes segfault
    #char  my_name[9], next_byte;
    #long  location;

    if DEBUG:
        print "op4_scan_b file=[%s]\n" %(filename)
    n_mat = 0;

    fp = open(filename, "rb");
    if DEBUG:
        print "op4_scan_b opened ok\n";

    while (not feof(fp)):
        location  = fp.tell();

        if DEBUG: printf("op4_scan_b location A =%ld\n", location);
        # Determining EOF is tricky because positioning the file pointer
        #  via Fortran record lengths ends up putting the FP at the very
        #  last byte in the file... without triggering EOF.  Have to
        #  explicitly look ahead one byte to get feof() to register .

        next_byte = fgetc(fp);
        if(feof(fp)): # ok I'm really done with this file
            break;
        else          # oops, am not at EOF, put the byte back
            ungetc(next_byte, fp);

        op4_is_mat_header_b(fp              ,
                            endian          ,
                           &record_length   ,
                           &is_header       ,
                            my_name         ,
                           &storage[n_mat] ,
                           &nRow[n_mat]    ,
                           &nCol[n_mat]    ,
                           &nStr[n_mat]    ,
                           &nNnz[n_mat]    ,
                           &type[n_mat]    ,
                           &form[n_mat]    ,
                           &offset[n_mat]);
        if DEBUG: print "op4_scan_b got RL=%d\n" %(record_length)
        fseek(fp, record_length + 2*BYTES_PER_WORD, SEEK_CUR) # skip header

        if DEBUG: print "op4_scan_b loc=%ld is_header=%d\n" %(fp.tell(), is_header)
        if is_header:
            strncpy(name[n_mat], my_name, 9);
            if storage[n_mat]: # sparse 1 or 2; need to count strings
                result = op4_count_str_b(fp, endian, type[n_mat],
                                         storage[n_mat],
                                         nRow[n_mat], nCol[n_mat],
                                        &nStr[n_mat], &nNnz[n_mat]);
                if not result:
                    fp.close()
                    n_mat = 0
                    return result
            ++(n_mat);

        if DEBUG:
            location  = fp.tell();
            print "op4_scan_b location B =%ld nStr=%d\n" %(location, nStr[n_mat-1]);

    fseek(fp, 0, SEEK_END);     # Set the file pointer to the end of file.
    offset[n_mat] = fp.tell();  # If this isn't done ftell() returns bogus.
    if DEBUG:
        print "op4_scan_b end byte position=%ld\n" %(offset[n_mat])
        fp.close()

    if not storage[n_mat-1]:
        # last matrix stored in dense scheme
        nNnz[n_mat-1] = nRow[n_mat-1]nCol[n_mat-1];
    return 1;



def op4_is_mat_header_b(fp,            # {{{1
                         endian    ,    # in  0=native   1=flipped
                         record_length, # out
                         is_header,   # out 1=is matrix hdr; 0=isn't
                         my_name ,    # out matrix name
                         storage   ,  # out 0=dn; 1=sp1; 2=sp2
                         nRow      ,  # out number of rows
                         nCol      ,  # out number of columns
                         nStr      ,  # out number of strings
                         nNnz      ,  # out number of nonzero terms
                         type      ,  # out 1=RS; 2=RD; 3=CS; 4=CD
                         form      ,  # out matrix form 6=symm, etc
                         offset):     # out byte offset to matrix
    """
    If this Fortran record contains a matrix header, reads the appropriate
    entries.  The file pointer is rewound to where it was when
    this routine was called.
    """
    DEBUG = 0;
    #int   is_ascii, nCol_, nRow_, Form_, Type_, rec_len,
    #      column_id, start_row, n_words, header[7];
    #long  location;
    #char name_ptr;

    is_header = 0;

    location   = fp.tell();
    fread(record_length, BYTES_PER_WORD, 1, fp);
    if DEBUG:
        print "op4_is_mat_header_b RL=%d location=%ld\n" %(record_length, location)
    if endian:
        record_length = flip_bytes_int(record_length)
if DEBUG:
    print "op4_is_mat_header_b RL=%d location=%ld\n" %(record_length, location)

    if(record_length == 24):
        # candidate for a matrix header
        fread(header, BYTES_PER_WORD, 7, fp); # + trailing rec marker
        # words header[4],header[5] make up the ASCII name--is it valid?
        name_ptr = header[4];
        is_ascii = 1;
        if(op4_valid_name(name_ptr)):
            strncpy(my_name, name_ptr,9)
        else:
            strncpy(my_name, "unnamed",9)

        nCol_  = header[0];
        nRow_  = header[1];
        Form_  = header[2];
        Type_  = header[3];
        if(endian):
            nCol_ = flip_bytes_int(nCol_)
            nRow_ = flip_bytes_int(nRow_)
            Form_ = flip_bytes_int(Form_)
            Type_ = flip_bytes_int(Type_)

        if DEBUG:
            print "op4_scan_b name=[%s] " %(my_name)
            print " is_ascii=%d nCols=%d nRows=%d Form=%d Type=%d\n" %(is_ascii,nCol_,nRow_,Form_,Type_)

        if(is_ascii and nCol_ > 0 and nRow_ and Form_ < 50 and Type_ < 10):
            # we have a winner; this is a matrix header
            is_header = 1
            nRow      = fabs(nRow_)
            nCol      = nCol_
            offset    = location
            type      = Type_
            form      = Form_
            nStr      = 0
            nNnz      = 0

            fread(&rec_len  , BYTES_PER_WORD, 1, fp)
            fread(&column_id, BYTES_PER_WORD, 1, fp)
            fread(&start_row, BYTES_PER_WORD, 1, fp)
            fread(&n_words  , BYTES_PER_WORD, 1, fp)
            if(endian):
                rec_len   = flip_bytes_int(rec_len)
                column_id = flip_bytes_int(column_id)
                start_row = flip_bytes_int(start_row)
                n_words   = flip_bytes_int(n_words)

            if DEBUG:
                print "op4_scan_b rl=%2d c=%2d r=%2d nw=%2d\n" %(rec_len, column_id, start_row, n_words)
            if(nRow_ < 0):
                storage = 2;   # sparse 2
            elif(not start_row or
                       ((start_row == 1) and (n_words == 1))): # null
                storage = 1;   # sparse 1
            else:              # dense
                storage = 0;

    fseek(fp, location, SEEK_SET);


def op4_count_str_b(FILE fp     ,   # in  {{{1
                     int   endian ,  # in
                     int   nType  ,  # in
                     int   storage,  # in
                     int   nRows  ,  # in
                     int   nCols  ,  # in
                     int  nStr   ,   # out
                     int  nNnz):     # out
    DEBUG = 0
    #int col, row, n_words, length_in_words, record_length, encoded, word_count, end_of_data;

    #define ErrStrSize 200
    #char error_msg[ErrStrSize - 1];
    nStr = 0
    col  = 1

    while (col <= nCols):

        if DEBUG:
            print "op4_count_str_b top of loop loc=%ld\n" %(fp.tell())
        fread(&record_length, BYTES_PER_WORD, 1, fp);
        fread(&col,           BYTES_PER_WORD, 1, fp);
        fread(&end_of_data,   BYTES_PER_WORD, 1, fp);
        fread(&n_words,       BYTES_PER_WORD, 1, fp);
        if endian:
            record_length = flip_bytes_int(record_length)
            col           = flip_bytes_int(col)
            end_of_data   = flip_bytes_int(end_of_data)
            n_words       = flip_bytes_int(n_words)

        if DEBUG:
            printf("record_length  =%d\n", record_length);
            printf("col            =%d\n", col          );
            printf("end_of_data    =%d\n", end_of_data  );
            printf("n_words        =%d\n", n_words      );

        # make sure values aren't bogus {{{2
        if((end_of_data != 0) and (end_of_data != 1)):
            nNnz = 0;
            snprintf(error_msg, ErrStrSize,
                    "end_of_data=%d out of range, byte offset %ld\n",
                     end_of_data, fp.tell());
            raise Exception(error_msg);

        if((col < 1) or (col > (nCols + 2))):
            nNnz = 0;
            snprintf(error_msg, ErrStrSize,
                    "col=%d out of range, byte offset %ld\n",
                     col, fp.tell());
            raise Exception(error_msg);

        if((record_length < 4*BYTES_PER_WORD) or
            (record_length > BYTES_PER_WORD *
                                (5 + op4_words_per_term[nType]nRows))):
                              #  5 = 3 column header ints +
                              #      2 string header ints
            nNnz = 0;
            snprintf(error_msg, ErrStrSize,
                    "RL=%d out of range (should be < %d), byte offset %ld\n",
                     record_length,
                     BYTES_PER_WORD*(5 + op4_words_per_term[nType]nRows),
                     fp.tell());

            raise Exception(error_msg);
        # 2}}}

        if DEBUG:
            printf("op4_count_str_b A RL=%d col=%d enddata=%d n_words=%d  loc=%ld\n" %(record_length, col, end_of_data, n_words, fp.tell());

        if end_of_data:
            fseek(fp, record_length - 2*BYTES_PER_WORD, SEEK_CUR);
            break;

        word_count = 0;
        while (word_count < n_words):
            if(storage == 1):
                fread(&encoded, BYTES_PER_WORD, 1, fp);
                word_count += 1;
                if endian:
                    encoded = flip_bytes_int(encoded);

                length_in_words = (encoded/65536) -1;
                if DEBUG:
                    print "op4_count_str_b B sp1 length_words=%d encoded=%d loc=%ld\n" %(length_in_words, encoded, fp.tell())
            else:  # storage == 2
                fread(&length_in_words, BYTES_PER_WORD, 1, fp);
                fread(&row,             BYTES_PER_WORD, 1, fp); # unused
                word_count += 2;
                if endian:
                    length_in_words = flip_bytes_int(length_in_words);
                
                --length_in_words;
            if DEBUG:
                printf("op4_count_str_b C sp2 length_words=%d row=%d  loc=%ld\n" %(length_in_words, row, fp.tell());
            
            if(length_in_words < 0):
                snprintf(error_msg, ErrStrSize,
                        "length_in_words=%d out of range, byte offset %ld\n",
                         length_in_words, fp.tell());
                raise Exception(error_msg);
            nNnz += length_in_words;
            ++(nStr);

            # skip numeric data
            fseek(fp, length_in_words*BYTES_PER_WORD, SEEK_CUR);
            word_count += length_in_words;

        if DEBUG:
            print "op4_count_str_b D word_count=%d  loc=%ld\n" %(word_count, fp.tell());
        
        fread(&record_length, BYTES_PER_WORD, 1, fp);
        if endian:
            record_length = flip_bytes_int(record_length);
        if DEBUG:
            print "op4_count_str_b E end RL=%d  loc=%ld\n" %(record_length, fp.tell())

    nNnz /= op4_words_per_term[nType];
    return 1
# 1}}}

def op4_line_type(line):  # in {{{1
    #int length, type;

    DEBUG = 0
    length = len(line)
    if(length > 2 and line[2] == '.'):  # most common case:  numeric data
        type = OP4_TEXT_NUMERIC;
    elif(length >= 34):  # header string must be at least 34 chars
        type = OP4_TEXT_HEADER;
    elif(length >= 24 and line[23] >= '0' and line[23] <= '9'):
        type = OP4_TEXT_NEW_COLUMN;
    elif(length >= 16 and line[15] >= '0' and line[15] <= '9'):
        type = OP4_TEXT_STRING_2;
    elif(length >=  8 and line[ 7] >= '0' and line[ 7] <= '9'):
        type = OP4_TEXT_STRING_1;
    else:
        if DEBUG:
            print "op4_line_type Error on [%s]\n" %(line)
        type = OP4_TEXT_ERROR
        if DEBUG:
            print "T=%d:[%s]\n", %(type,line)
    return type;
# 1}}}

def op4_read_col_t(fp         ,  # {{{1
                   c_in       ,  # in  requested column to read
                   nRow       ,  # in  # rows    in matrix
                   nCol       ,  # in  # columns in matrix
                   fmt_str    ,  # in  eg "%23le%23le%23le"
                   col_width  ,  # in  # characters in format str
                   storage    ,  # in  0=dense  1,2=sparse
                   complx     ,  # in  0=real   1=complex
                   n_str      ,  # in/out idx str_data[] (s_o)=1,2
                   S          ,  # out string data   (s_o)=1,2
                   N_index    ,  # in/out idx N[]    (s_o)=1,2
                   N             # out numeric data
    DEBUG = 0
    #double      x[OP4_MAX_TEXT_COLS];
    #int         nRead, j, nNum_cols, still_reading, location, col, row,
    #            nNum, got_this_col, n_nnz, mixed, nwords, start_row,
    #            line_type, NPT,
    #            n_read_this_str = 0;
    #char  line[OP4_TXT_LINE_SIZE + 1];

    still_reading = 1;
    got_this_col  = 0;
    n_nnz         = 0;  # count of numeric values loaded
    NPT           = Nm_per_term[complx]; # numbers per term; real=1 complex=2
    if not storage:
        # zero out dense column
        for j in range(NPTnRow):
            N[j] = 0.0;

    while (still_reading):
        location = fp.tell();
        fgets(line, OP4_TXT_LINE_SIZE, fp);
        line_type = op4_line_type(line);
        if DEBUG:
            print "op4_read_col_t line_type=%d\n" %(line_type)
        if(strlen(line) >= 8):
            switch (line_type)
                case OP4_TEXT_NUMERIC   : # {{{2
                    nNum_cols = (int) (len(line) / col_width)
                    nRead = sscanf(line, fmt_str, &x[0], &x[1], &x[2],
                                                  &x[3], &x[4], &x[5],
                                                  &x[6], &x[7], &x[8]);
                    if storage:  # sparse
                        if DEBUG:
                            print "op4_read_col_t inserting %d terms at row=%d\n" %(nRead, S[n_str-1].start_row)
                        n_read_this_str += nRead;
                        S[n_str-1].len = n_read_this_str/NPT;
                        for j in range(nRead):
                            N[(N_index)++] = x[j];
                    else:       # dense
                        if DEBUG: print "op4_read_col_t inserting %d terms at row=%d\n" %(nRead, row);
                        for j in range(nRead):
                            if DEBUG:
                                print " N[row=%d + N_index=%d]=%e\n" %(row, (N_index), x[j])
                            N[row + j] = x[j];
                        row += nRead;

                    n_nnz += nRead;
                    if DEBUG:
                        print "op4_read_col_t nR=%d n_nnz now=%d\n" %(nRead, n_nnz);
                case OP4_TEXT_HEADER    : # {{{2
                    raise Exception("op4_read_col_t new matrix header while reading column");
                case OP4_TEXT_NEW_COLUMN: # {{{2
                    if DEBUG:
                        print "op4_read_col_t start new column (got_this_col=%d) loc=%ld\n" %(got_this_col, fp.tell());
                    if(got_this_col):
                        # have read next column's header; back up & exit
                        fseek(fp, location, SEEK_SET);
                        still_reading = 0;
                    else:
                        nRead = sscanf(line, "%8d%8d%8d", &col, &row, &nNum);
                        if DEBUG:
                            print " requested col=%d  this col=%d  this row=%d\n", c_in, col, row)
                        if(c_in != col):
                            # this is not the column of interest
                            fseek(fp, location, SEEK_SET);
                            return n_nnz;
                        --col; --row; # internal indexing is 0 based
                        if not storage: # dense
                            row *= NPT
                        if DEBUG:
                            print "\n new column %2d row %2d complx=%d\n" %(col, row, complx);
                        got_this_col = 1;
                        if((col > nCol) or (nRead != 3)):
                            # have reached the spurious end of matrix entry
                            # or have encountered bad data
                            return n_nnz;

                case OP4_TEXT_STRING_1  :
                case OP4_TEXT_STRING_2  : # {{{2
                    n_read_this_str = 0;
                    if(line_type == OP4_TEXT_STRING_1):
                        sscanf(line, "%8d", &mixed);
                        nwords    = (mixed/65536) - 1;
                        start_row = mixed - 65536 * (nwords + 1);
                    else: # type 2
                        sscanf(line, "%8d%8d", &nwords, &start_row);
                        --nwords;

                    S[n_str].start_row = start_row - 1;
                    S[n_str].len       = 0;
                    S[n_str].N_idx     = N_index;
                    ++(n_str);
                case OP4_TEXT_ERROR: # {{{2
                    raise Exception("op4_read_col_t error reading a column ");

    return n_nnz;
# 1}}}

def op4_read_col_b(fp         ,  # {{{1
                   endian     ,  # in  0=native   1=flipped
                   c_in       ,  # in  requested column to read
                   nRow       ,  # in  # rows    in matrix
                   nCol       ,  # in  # columns in matrix
                   nType      ,  # in  1=RS 2=RD 3=CS 4=CD
                   storage    ,  # in  0=dn; 1=sp1; 2=sp2
                   n_str      ,  # in/out idx str_data[] (s_o)=1,2
                   S          ,  # out string data   (s_o)=1,2
                   N_index    ,  # in/out idx N[]    (s_o)=1,2
                   N             # out numeric data
                   ):
    DEBUG = 0
    #int    complx, NPT, WPN, BPN, i, record_length, col, start_row,
    #       n_numbers, got_this_col, n_nnz, is_header, unused, mixed,
    #       n_w_this_col, n_w_this_str;
    #long   l_unused;
    #char   unused_name[9];
    #float  xs;

    if DEBUG:
        print "\nop4_read_col_b 1 location=%ld  endian=%d want col %d\n" %(fp.tell(), endian, c_in);

    complx        = nType > 2 ? 1 : 0;
    NPT           = Nm_per_term[complx];  # numbers per term; real=1 complex=2
    WPN           = nType % 2 ? 1 : 2;    # words per number; float=1 double=2
    BPN           = BYTES_PER_WORD * WPN; # bytes per number
    got_this_col  = 0;
    n_nnz         = 0;  # count of numeric values loaded
    op4_is_mat_header_b(fp              , # in
                        endian          , # in
                       &record_length   , # out
                       &is_header       , # out
                        unused_name     , # out
                       &unused          , # out
                       &unused          , # out
                       &unused          , # out
                       &unused          , # out
                       &unused          , # out
                       &unused          , # out
                       &unused          , # out
                       &l_unused);        # out
    if DEBUG:
        print "op4_read_col_b 2 RL=%d location=%ld\n" %(record_length, fp.tell())
    if(is_header): # yes?  skip it
        fseek(fp, record_length+2*BYTES_PER_WORD, SEEK_CUR);
    if DEBUG:
        print "op4_read_col_b 3 RL=%d location=%ld\n" %(record_length, fp.tell())

    if not storage:
        # zero out dense column
        for i in range(NPTnRow):
            N[i] = 0.0;

    fread(&record_length, BYTES_PER_WORD, 1, fp);
    fread(&col,           BYTES_PER_WORD, 1, fp);
    fread(&start_row,     BYTES_PER_WORD, 1, fp);
    fread(&n_w_this_col,  BYTES_PER_WORD, 1, fp);
    if endian:
        record_length = flip_bytes_int(record_length);
        col           = flip_bytes_int(col);
        start_row     = flip_bytes_int(start_row);
        n_w_this_col  = flip_bytes_int(n_w_this_col);

    if DEBUG:
        print "op4_read_col_b: requested col=%d  RL=%d got col=%d row=%d nW=%d\n" %(c_in, record_length, col, start_row, n_w_this_col)
    n_numbers = n_w_this_col/WPN;

    if(c_in > col): # requested column ahead of me; skip record
        fseek(fp, record_length-2*BYTES_PER_WORD, SEEK_CUR);
        if DEBUG:
            print "op4_read_col_b: column mismatch, advance to location %ld & exit\n" %(fp.tell());
        return 0;

    elif(c_in < col): # requested column behind me; back up
        fseek(fp,-4*BYTES_PER_WORD, SEEK_CUR);
        if DEBUG:
            print "op4_read_col_b: column mismatch, back up to location %ld & exit\n" %(fp.tell())
        return 0;

    if DEBUG:
        print "op4_read_col_b: going to read col=%d\n" %(col)

    if not storage:             # dense; read entire column
        if(WPN == 2):
            # double precision:  can do bulk read
            if DEBUG: print "op4_read_col_b: double prec\n"
            fread(&N[NPT*(start_row-1)], BPN, n_numbers, fp);
            if endian: # opposite endian, need to flip terms
                for i in range(n_numbers):
                    N[NPT*(start_row-1) + i] = flip_bytes_double(N[NPT*(start_row-1) + i]);

        elif not endian:
            # single precision native endian
            if DEBUG: print "op4_read_col_b: single prec native endian\n"
            for i in range(n_numbers):
                fread(&xs, BPN, 1, fp);
                N[NPT*(start_row-1) + i] = xs;
        else:
            # single precision opposite endian
            if DEBUG: print "op4_read_col_b: single prec opposite endian\n"
            for i in range(n_numbers):
                fread(&xs, BPN, 1, fp);
                N[NPT*(start_row-1) + i] = flip_bytes_float(xs);

        if DEBUG:
            for i in range(n_numbers):
                print "op4_read_col_b: N[%d]=%le\n" %(NPT*(start_row-1) + i,N[NPT*(start_row-1) + i])
        fseek(fp, BYTES_PER_WORD, SEEK_CUR);  # skip trailing RL

        if DEBUG:
            print "op4_read_col_b: finished numeric data location=%ld\n" %(fp.tell());

    else:  # sparse; loop over all strings
        record_length -= 2*BYTES_PER_WORD; # don't count start, end rec mark
        while (record_length >= 8):
            if(storage == 1):  # sparse 1
                fread(&mixed, BYTES_PER_WORD, 1, fp);
                if endian:
                    mixed = flip_bytes_int(mixed);

                n_w_this_str  = (mixed/65536) - 1;
                start_row     = mixed - 65536 * (n_w_this_str + 1);
                record_length -= 1*BYTES_PER_WORD;
            else:             # sparse 2
                fread(&n_w_this_str, BYTES_PER_WORD, 1, fp);
                fread(&start_row,    BYTES_PER_WORD, 1, fp);
                record_length -= 2*BYTES_PER_WORD;
                if endian:
                    n_w_this_str = flip_bytes_int(n_w_this_str);
                    start_row    = flip_bytes_int(start_row);
                --n_w_this_str;
            n_numbers = n_w_this_str/WPN;
            if DEBUG:
                print "\nop4_read_col_b top of str: start_row=%d n_w_this_str=%d n_numbers=%d N_index=%d\n" %(start_row, n_w_this_str, n_numbers, N_index)

            S[n_str].start_row = start_row - 1;
            S[n_str].len       = n_numbers/(complx + 1);
            S[n_str].N_idx     = N_index;
            ++(n_str);
            if(WPN == 2):
                # double precision:  can do bulk read
                if DEBUG: print "op4_read_col_b sp: double prec\n"

                fread(&N[N_index], BPN, n_numbers, fp);
                if endian: # opposite endian, need to flip terms
                    for i in range(n_numbers):
                        N[N_index + i] = flip_bytes_double(N[N_index + i]);
            elif not endian:
                # single precision native endian
                if DEBUG: print "op4_read_col_b sp: single prec native endian\n")
                for i in range(n_numbers):
                    fread(&xs, BPN, 1, fp);
                    N[N_index + i] = xs;
            else:
                # single precision opposite endian
                if DEBUG: print "op4_read_col_b sp: single prec opposite endian\n"
                for i in range(n_numbers):
                    fread(&xs, BPN, 1, fp);
                    N[N_index + i] = flip_bytes_float(xs);
            record_length -= BYTES_PER_WORD*n_w_this_str;
            if DEBUG:
                print "op4_read_col_b sp end string  RL=%d loc=%ld\n" %(record_length, fp.tell());
                for i in range(n_numbers):
                    print "op4_read_col_b N[%3d] = %e\n" %(N_index+i, N[N_index + i])
            N_index += n_numbers
            n_nnz   += n_numbers
        fread(&record_length, BYTES_PER_WORD, 1, fp)
    return n_nnz;

FILE* op4_open_r(const char filename, long offset):  # in
    FILE fp;
    # Open an op4 file for reading, optionally positioning the file
    # handle to a location past the start of the file.
    #

    fp = fopen(filename, "r");
    if(fp == NULL)
        return NULL;
    if(offset)
        fseek(fp, offset, SEEK_SET);
    return fp;

def int flip_bytes_int(int x):
    #int   y;
    #char *c_x, *c_y;

    c_x = (char *) &x;
    c_y = (char *) &y;

    c_y[0] = c_x[3];
    c_y[1] = c_x[2];
    c_y[2] = c_x[1];
    c_y[3] = c_x[0];

    return y;

def float flip_bytes_float(float x):
    """
    Looks identical to flip_bytes_int() and it practice it might
    be possible to use this function and flip_bytes_int()
    interchangably.  The combination of gcc and x86 hardware may
    cause problems though since floats might be assigned extra
    bytes internally.
    """
    float   y;
    char    *c_x, *c_y;

    c_x = (char *) &x;
    c_y = (char *) &y;

    c_y[0] = c_x[3];
    c_y[1] = c_x[2];
    c_y[2] = c_x[1];
    c_y[3] = c_x[0];

    return y;

def double flip_bytes_double(double x):
    double   y;
    char    *c_x, *c_y;

    c_x = (char *) &x;
    c_y = (char *) &y;

    c_y[0] = c_x[7];
    c_y[1] = c_x[6];
    c_y[2] = c_x[5];
    c_y[3] = c_x[4];
    c_y[4] = c_x[3];
    c_y[5] = c_x[2];
    c_y[6] = c_x[1];
    c_y[7] = c_x[0];

    return y;

def op4_valid_name(name):
     """ Returns:
         1 if the name passed in is a valid Nastran op4 matrix name
        -1 if the name is null
         0 if the name is not valid
     @todo this check is too conservative; Nastran allows names like 'a+'
     """
    DEBUG = 0

    #int i, len, fixed_len;
    if DEBUG:
        print "op4_valid_name with [%s]\n" %(name)

    length = len(name)

    # strip trailing blanks
    fixed_len = len;
    name = name.strip()
    len = len(name)

    if not length:                    # null string?
        return -1;
    elif(length > 8):               # too long?
        return  0;
    elif(not isalpha((int) name[0])):  # first char must be [a-zA-Z]
        return  0;

    for i in range(len):
        # this test is too conservative; Nastran allows names like 'a+'
        if(not (ISUNERDSCORE( name[i]) or
              isalpha((int) name[i]) or
              isdigit((int) name[i]))):
            return 0;

    return 1;

# words



def mexFunction(nlhs, plhs[], int nrhs, prhs[]):
    DEBUG = 0
    filename[80],
          line[OP4_TXT_LINE_SIZE+1],
          name[Max_Matrices_Per_OP4][9],
          fmt_str[OP4_FMT_LINE_SIZE],
          fmt_one[OP4_FMT_STR_SIZE];
    int   c, i, j, s_ptr, n_ptr, n_nnz, n_mat_in_file, complx, rows, NPT,
          Nmat, Nskip, nNum_cols, filetype,
          col_width = 0,
          load_all  = 0,
          unused    = 0,
          storage[Max_Matrices_Per_OP4],
          nRow[Max_Matrices_Per_OP4],
          nCol[Max_Matrices_Per_OP4],
          nStr[Max_Matrices_Per_OP4],
          nNnz[Max_Matrices_Per_OP4],
          type[Max_Matrices_Per_OP4],
          form[Max_Matrices_Per_OP4],
          digits[Max_Matrices_Per_OP4];
    long  offset[Max_Matrices_Per_OP4];

    #double  nskip_ptr;
    double  Mr[Max_Matrices_Per_OP4], Mi[Max_Matrices_Per_OP4];
    int     ir[Max_Matrices_Per_OP4], jc[Max_Matrices_Per_OP4];
    #int     r, s, ir_index, n_str_one_col,
    #        ptr_H, ptr_S_start, ptr_N_start, ptr_S, ptr_N, size
    print_header = 1;  # 1=print matrix info after loading; 0=don't

    SparseMatrix m;
    D = 0.0;

    if((nrhs < 0) or (nrhs > 2)):
        raise Exception("Can only have one or two RHS arguments.\n")

    if(not mxIsChar(prhs[0])):
        raise Exception("First RHS argument must be a text filename.\n")

    mxGetString(prhs[0], filename, 80)
    Nmat = nlhs
    if((Nmat < 1) or (Nmat >= Max_Matrices_Per_OP4) ):
        raise Exception("Number of LHS arguments must be >= 1 and <= 50.\n")

    if(nrhs > 1):     # get # of matrices to skip over
        nskip_ptr = mxGetPr(prhs[1])
        Nskip     = (int) nskip_ptr
        if(Nskip < 0):
            raise Exception("Second RHS argument must be an integer >= 0.\n")
    else:
        Nskip = 0


    if DEBUG:
        print "loadop4 nskip=%d  nmat=%d  name=[%s]\n" %(Nskip, Nmat, filename);

    # Scan the file for number of matrices and headers for each matrix.
    if(not op4_scan(filename,      # in
                 &n_mat_in_file, # out number of matrices
                  name,          # out matrix names
                  storage,       # out 0=dn; 1=sp1; 2=sp2
                  nRow,          # out number of rows
                  nCol,          # out number of columns
                  nStr,          # out number of strings
                  nNnz,          # out number of nonzero terms
                  type,          # out 1=RS; 2=RD; 3=CS; 4=CD
                  form,          # out matrix form 6=symm, etc
                  digits,        # out size of mantissa
                  offset)):      # out byte offset to matrix
        snprintf(line, OP4_TXT_LINE_SIZE, " file error for '%s'", filename);
        raise Exception(line);

    if DEBUG:
        print "loadop4 n_mat_in_file = %d\n" %(n_mat_in_file)
    if(n_mat_in_file < 1 ):
        snprintf(line, OP4_TXT_LINE_SIZE, " no matrices in '%s'", filename);
        raise Exception(line);
        return

    if(Nmat == -1):  # wants every matrix in the file (beyond nSkip)
        load_all = 1
        Nmat     = n_mat_in_file - Nskip

    if((Nmat < 0) or (Nskip + Nmat) > n_mat_in_file):
        # Requested more matrices than exist in file.
        snprintf(line, OP4_TXT_LINE_SIZE, " only %d matrices in '%s'",n_mat_in_file, filename);

        raise Exception(line);

    filetype = op4_filetype(filename);

    fp = op4_open_r(filename, 0);
    if(fp == NULL):
        raise Exception("loadop4: unable to open op4 file at offset\n");

    for i in range(Nskip,Nskip+Nmat):
        if DEBUG:
            print "%2d %2d x %2d %2s %2s nStr=%4d nNnz=%4d f=%2d %-8s D=%2d O=%12ld L=%ld\n" %(i+1, nRow[i], nCol[i], op4_type_str[type[i]-1],
                                                                                               op4_store_str[storage[i]], nStr[i], nNnz[i], form[i], name[i],
                                                                                               digits[i], offset[i], offset[i+1] - offset[i])
        # allocate memory, place stack item to receive op4 matrix {{{2
        complx = type[i] > 2 ? 1 : 0;
        NPT = Nm_per_term[complx];

        if storage[i]: # sparse
            # matlab sparse {{{3

            # Declare a tops native sparse matrix having only 1 column.
            # Assume worst case values for # strings, # non-zeros.  Have
            # to declare a double to force alignment to multiple of 8.
            # This code duplicates code in the following functions:
            #  - src/sparce.c:sp_data_offsets()
            #  - src/sparce.c:sp_set_header()
            #  - src/sparce.c:sparse_overlay()
            n_str_one_col = nRow[i]/2 + 1;
            ptr_H         =               SPARSE_MAGIC_SIZE;
            ptr_S_start   = ptr_H       + SPARSE_HDR_SIZE;
            ptr_N_start   = ptr_S_start +     (1 + 1);  # nCOl = 1
            ptr_S         = ptr_N_start +     (1 + 1);  # nCOl = 1
            ptr_N         = ptr_S       + n_str_one_col*sizeof(str_t)/
                                                        sizeof(int);

            if(ptr_N % 2):
                ++(ptr_N)

            size = (ptr_N)*sizeof(int) + NPTnRow[i]*sizeof(double)

            if((D = (double *) malloc(size)) == NULL):
                raise Exception("loadop4:  insufficient memory for matrix column");

            m.data        = (char   *)   D;
            m.magic       = (int *) m.data;
            m.H           =            &m.magic[SPARSE_MAGIC_SIZE];
            m.S_start     =            &m.magic[ptr_S_start];
            m.N_start     =            &m.magic[ptr_N_start];
            m.S           = (str_t  *) &m.magic[ptr_S      ];
            m.N           = (double *) &m.magic[ptr_N      ];

            m.H[ROWS]     = nRow[i];
            m.H[COLS]     = 1;
            m.H[COMPLX]   = complx;
            m.H[n_STR]    = nRow[i]/2 + 1;
            m.H[n_NONZ]   = nRow[i]*NPT;
            m.H[DATA_SIZE]= size;

            # then create the matlab variables
            if complx:           # sparse complex
                plhs[i-Nskip] = mxCreateSparse(nRow[i], nCol[i],
                                               nNnz[i], mxCOMPLEX);
                Mr[i-Nskip]   = mxGetPr(plhs[i-Nskip]);
                Mi[i-Nskip]   = mxGetPi(plhs[i-Nskip]);
            else:                # sparse real
                plhs[i-Nskip] = mxCreateSparse(nRow[i], nCol[i],
                                               nNnz[i], mxREAL);
                Mr[i-Nskip]   = mxGetPr(plhs[i-Nskip]);
            ir[i-Nskip] = mxGetIr(plhs[i-Nskip]);
            jc[i-Nskip] = mxGetJc(plhs[i-Nskip]);

            # 3}}}

        else:          # dense
            if complx:         # dense complex
                plhs[i-Nskip] = mxCreateDoubleMatrix(nRow[i], nCol[i],
                                                     mxCOMPLEX);
                Mr[i-Nskip]   = mxGetPr(plhs[i-Nskip]);
                Mi[i-Nskip]   = mxGetPi(plhs[i-Nskip]);
                for j in range(nRow[i]*nCol[i]):
                    Mi[i-Nskip][j] = 0.0;
            else:              # dense real
                plhs[i-Nskip] = mxCreateDoubleMatrix(nRow[i], nCol[i],
                                                     mxREAL);
                Mr[i-Nskip]   = mxGetPr(plhs[i-Nskip]);
            for j in range(nRow[i]*nCol[i]):
                Mr[i-Nskip][j] = 0.0;

        # 2}}}

        fseek(fp, offset[i], SEEK_SET);  # jump to the ith matrix

        if(filetype == 1):  # text
            fgets(line, OP4_TXT_LINE_SIZE, fp);  # skip the header line
            # create the scanf format string, eg
            col_width =  digits[i] + 7;
            strncpy( fmt_str, "", OP4_FMT_STR_SIZE);
            snprintf(fmt_one, OP4_FMT_STR_SIZE, "%%%dle", col_width);
            nNum_cols = (int) (80 / col_width);
            for c in range(nNum_cols):
                strncat(fmt_str, fmt_one, OP4_FMT_STR_SIZE);
            if DEBUG:
                print "one=[%s] str=[%s]\n" %(fmt_one, fmt_str)

        if storage[i]:   # sparse
            jc[i-Nskip][0] = 0;
            ir_index       = 0;
        else:            # dense
            if((D = (double *) malloc(nRow[i]*NPT*sizeof(double))) == NULL):
                raise Exception("loadop4:  insufficient memory for matrix column");

        for c in range(1,nCols[i]):  # op4 columns start at 1
            if DEBUG:
                print "loadop4->column %2d n_ptr=%d\n" %(c, n_ptr)

            # read column c {{{2
            if storage[i]:     # sparse
                s_ptr            = 0;
                n_ptr            = 0;
                m.S_start[0]     = 0;
                m.N_start[0]     = 0;
                m.S[0].start_row = 0;
                m.S[0].len       = 0;
                m.S[0].N_idx     = 0;

            if(filetype == 1):  # text
                if storage[i]:     # load as sparse matrix
                    n_nnz = op4_read_col_t(fp, c, nRow[i], nCol[i], fmt_str, col_width,
                                   storage[i], # in 0=dn  1,2=sp1,2  3=ccr
                                   complx    , # in  0=real   1=complex
                                  &s_ptr     , # in/out index m.S (if sp1,2)
                                   m.S       , # out string data (if sp1,2)
                                  &n_ptr     , # in/out index m.N
                                   m.N);       # out numeric data

                else:                   # load as dense matrix
                    n_nnz =
                    op4_read_col_t(fp, c, nRow[i], nCol[i], fmt_str, col_width,
                                   storage[i], # in 0=dn  1,2=sp1,2  3=ccr
                                   complx    , # in  0=real   1=complex
                                   unused    , # in/out index m.S (if sp1,2)
                         (str_t *) unused    , # out string data  (if sp1,2)
                         (int   *) unused    , # in/out index m.N (sp 1,2)
                                   D); # out numeric data
            else:              # binary
                if storage[i]:     # load as sparse matrix
                    n_nnz = op4_read_col_b(fp, filetype - 2,
                                   c, nRow[i], nCol[i], type[i],
                                   storage[i], # in 0=dn  1,2=sp1,2
                                  &s_ptr     , # in/out index m.S (if sp1,2)
                                   m.S       , # out string data  (if sp1,2)
                                  &n_ptr     , # in/out index m.N
                                   m.N);       # out numeric data

                else:                   # load as dense matrix
                    n_nnz =
                    op4_read_col_b(fp, filetype - 2,
                                   c, nRow[i], nCol[i], type[i],
                                   storage[i], # in 0=dn  1,2=sp1,2
                                   unused    , # in/out index m.S (if sp1,2)
                         (str_t *) unused    , # out string data (if sp1,2)
                         (int   *) unused    , # in/out index m.N (sp 1,2)
                                   D
                                  ); # out numeric data
            # # 2}}}


            # copy column c from internal to matlab variable
            if storage[i]:   # sparse
                for s in range(s_ptr):
                    for j in range(m.S[s].len):
                        ir[i-Nskip][ir_index] = m.S[s].start_row + j;
                        Mr[i-Nskip][ir_index] = m.N[m.S[s].N_idx + j*NPT];
                        if complx:
                            Mi[i-Nskip][ir_index] = m.N[m.S[s].N_idx + j*NPT+1];
                        ++ir_index;
                jc[i-Nskip][c] = ir_index;

            else: # dense
                for r in range(nRow[i]):
                    Mr[i-Nskip][(c-1)nRow[i] + r] = D[r*NPT];
                    if complx:
                        Mi[i-Nskip][(c-1)nRow[i] + r] = D[r*NPT + 1];

            if DEBUG:
                print "loadop4<-column %2d read %d values\n\n" %(c, n_nnz)
        # end loop over columns of matrix i

        if(not storage[i])     # dense matrix
            free(D);
        if(print_header):
            print "loadop4: %2d x %2d %2s %2s nStr=%4d nNnz=%4d f=%2d %-8s dg=%1d of %1d\n" %(nRow[i], nCol[i], op4_type_str[type[i]-1],op4_store_str[storage[i]], nStr[i], nNnz[i],form[i],name[i], digits[i], i+1, Nmat)
    # # end loop over matrices

    fp.close()
    if DEBUG:
        print "end of loadop4\n"
    return;
